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Abstract 

We present an analytical solution and numerical tests of the epidemic-type aftershock (ETAS) 
model for aftershocks, which describes foreshocks, aftershocks and mainshocks on the same 
footing. In this model, each earthquake of magnitude m triggers aftershocks with a rate 
proportional to 10 am . The occurrence rate of aftershocks triggered by a single mainshock 
decreases with the time from the mainshock according to the modified Omori law Kj (t + c) p 
with p — 1 + 9. Contrary to the usual definition, the ETAS model does not impose an aftershock 
to have a magnitude smaller than the mainshock. Starting with a mainshock at time t = that 
triggers aftershocks according to the local Omori law, that in turn trigger their own aftershocks 
and so on, we study the seismicity rate of the global aftershock sequence composed of all the 
secondary and subsequent aftershock sequences. The effective branching parameter n, defined as 
the mean aftershock number triggered per event, controls the transition between a sub-critical 
regime n < 1 to a super-critical regime n > 1. A characteristic time t* , function of all the ETAS 
parameters, marks the transition from the early time behavior to the large time behavior. In the 
sub-critical regime, we recover and document the crossover from an Omori exponent 1 — 9 for 
t < t* to 1 + 9 for t > t* found previously in [Sornette and Sornette, 1999] for a special case of 
the ETAS model. In the super-critical regime n > 1 and 9 > 0, we find a novel transition from 
an Omori decay law with exponent 1 — 9 for t < t* to an explosive exponential increase of the 
seismicity rate for t > t*. The case 9 < yields an infinite n- value. In this case, we find another 
characteristic time r controlling the crossover from an Omori law with exponent 1 — \9\ for t < r, 
similar to the local law, to an exponential increase at large times. These results can rationalize 
many of the stylized facts reported for aftershock and foreshock sequences, such as (i) the 
suggestion [Liu, 1984; Bowman, 1997] that a small p-value may be a precursor of a large 
earthquake, (ii) the relative seismic quiescence sometimes observed before large aftershocks, (iii) 
the positive correlation between b and p-values, (iv) the observation that great earthquakes are 
sometimes preceded by a decrease of 6-value and (v) the acceleration of the seismicity preceding 
great earthquakes. 
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Introduction 

It is well known that the seismicity rate increases 
after a large earthquake, for time period up to one 
hundred years [Utsu et al, 1995], and distances up 
to several hundred km [Tajima and Kanamori, 1985; 
Steeples and Steeples, 1996; Kagan and Jackson, 1998; 
Meltzner and Wald, 1999; Dreger and Savage, 1999]. 
The rate of the triggered events usually decays in time 
as the modified Omori law n(t) — K/(t + c) p , where 
the exponent p is found to vary between 0.3 and 2 
[Davis and Frohlich, 1991; Kisslinger and Jones, 1991; 
Guo and Ogata, 1995; Utsu et al., 1995] and is often 
close to 1 (see however [Kisslinger, 1993; Gross and 
Kisslinger, 1994] for alternative decay laws such as 
the stretched exponential). 

These triggered events are called aftershocks if 
their magnitude is smaller than the first event. How- 
ever, the definition of an aftershock contains unavoid- 
ably a degree of arbitrariness because the qualification 
of an earthquake as an aftershock requires the spec- 
ification of time and space windows. In this spirit, 
several alternative algorithms for the definition of af- 
tershocks have been proposed [Gardner and Knopoff, 
1974; Molchan and Dmitrieva, 1992] and there is no 
consensus. 

Aftershocks may result from several and not nec- 
essarily exclusive mechanisms (see [Harris, 2001] and 
references therein): pore-pressure changes due to 
pore-fluid flows coupled with stress variations, slow 
redistribution of stress by aseismic creep, rate-and- 
statc dependent friction within faults, coupling be- 
tween the viscoelastic lower crust and the brittle up- 
per crust, stress-assisted micro-crack corrosion [Ya- 
mashita and Knopoff, 1987; Lee and Sornette, 2000], 
slow tectonic driving of a hierarchical geometry with 
avalanche relaxation dynamics [Huang et. al, 1998], 
dynamical hierarchical models with heterogeneity, feed- 
backs and healing [Blanter et al, 1997], etc. 

Since the underlying physical processes are not 
fully understood, the qualifying time and space win- 
dows are more based on common sense than on hard 
science. Particularly, there is no agreement about the 
duration of the aftershock sequence and the maximum 
distance between aftershock and mainshock. If one 
event occurs with a magnitude larger than the first 
event, it becomes the new mainshock and all preced- 
ing events are retrospectively called foreshocks. Thus, 
there is no way to identify foreshocks from usual af- 
tershocks in real time. There is also no way to distin- 
guish aftershocks from individual earthquakes [Hough 



and Jones, 1997]. The aftershock magnitude distribu- 
tion follows the Gutenberg-Richter distribution with 
similar &-value as other earthquakes [Ranalli, 1969; 
Knopoff et al, 1982]. They have also similar rupture 
process. Moreover, an event can be both an after- 
shock of a preceding large event, and a mainshock of 
a following earthquake. For example, the M=6.5 Big 
Bear event is usually considered as an aftershock of 
the M=7.3 Landers event, and has clearly triggered 
its own aftershock sequence. One can trace the dif- 
ficulty of the problem from the long-range nature of 
the interactions between faults in space and time re- 
sulting in a complex self-organized crust. 

In view of the difficulties in classifying sometimes 
an earthquake as a foreshock, a mainshock or an af- 
tershock, it is natural to investigate models in which 
this distinction is removed and to study their possible 
observable consequences. In this spirit, the epidemic 
type aftershock (ETAS) model introduced by Kagan 
and Knopoff [1981, 1987] and Ogata [1988] provides a 
tool for understanding the temporal clustering of the 
seismic activity without distinguishing between after- 
shocks, foreshocks and mainshock events. The ETAS 
model is a generalization of the modified Omori law, 
which takes into account the secondary aftershocks 
sequences triggered by all events. In this model, 
all earthquakes are simultaneously mainshocks, after- 
shocks and possibly foreshocks. An observed "after- 
shock" sequence is in the ETAS model the result of 
the activity of all events triggering events triggering 
themselves other events, and so on, taken together. 
The ETAS model aims at modeling complex after- 
shocks sequences and global seismic activity. The seis- 
micity rate is given by the superposition of aftershock 
sequences of all events. Each earthquake of magni- 
tude m triggers aftershock with a rate proportional 
to 10 Qm with the same coefficient a for all earth- 
quakes. The occurrence rate of aftershocks decreases 
with the time from the mainshock according to the 
modified Omori law K/(t + cf . The background seis- 
micity rate is modeled by a stationary Poisson process 
with a constant occurrence rate \i. Contrary to the 
usual definition, the ETAS model does not impose 
an aftershock to have a magnitude smaller than the 
mainshock. This way, the same law describes both 
foreshocks, aftershocks and mainshocks. This model 
has been used to give short-term probabilistic fore- 
cast of seismic activity [Kagan and Knopoff, 1987; 
Kagan and Jackson, 2000; Console and Murru, 2001], 
and to describe the temporal and spatial clustering of 
seismic activity [Ogata, 1988, 1989, 1992, 1999, 2001; 
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Kagan, 1991; Felzer et al, 2001]. Allthough the el- 
ementary results on the stability of the process have 
been known for many years [Kagan, 1991], no attempt 
has been made to study this model analytically in or- 
der to characterize its different regimes and obtain a 
deeper understanding of the combined interplay be- 
tween the model parameters (b, a, p, K, c and n) on 
the seismic activity. We stress below the contrast be- 
tween previous works in the mathematical statistical 
literature and our results. 

It should be noted that the ETAS model suf- 
fers from an important defect: it is fundamentally a 
"branching" model [Harris, 1963; Vere- Jones, 1977], 
with no "loops" . What this means is that an event 
has a unique "mother-mainshock" and not several. 
In the real case, we can expect that some events may 
be triggered by the combined loading and action at 
distance in time and space of several previous earth- 
quakes. Hence, events should have several "mothers" 
in general. This neglecting of "loops" is known in sta- 
tistical physics as a "mean-field" approximation and 
allows us to simplify the analysis while still keeping 
the essential physics in a qualitative way, even if the 
details may not be precisely recovered quantitatively. 

Sornette and Sornette [1999] studied analytically a 
particular case of the ETAS model, in which the af- 
tershock number does not depend on the mainshock 
magnitude, i.e., for a = 0. Starting with one event at 
time t — and considering that each earthquake gen- 
erates an aftershock sequence with a "local" Omori 
exponent p = 1 + 9, where 9 is a positive constant, 
they studied the decay law of the "global" aftershock 
sequence, composed of all secondary aftershock se- 
quences. They found that the global aftershock rate 
decays according to an Omori law with an exponent 
p = 1 — 6, smaller than the local one, up to a char- 
acteristic time t* , and then recovers the local Omori 
exponent p = 1 + 9 for time larger than t* . 

Here, we generalize their analysis in the more gen- 
eral case a > of the ETAS model, which includes a 
realistic magnitude distribution. We study the decay 
law of the global aftershock sequence as a function 
of the model parameters (local Omori law parame- 
ters and magnitude distribution). In addition to giv- 
ing more complete analytical results, we present nu- 
merical simulations that test these predictions. We 
also generalize the investigation and analysis into the 
"super-critical" regime. Indeed, depending on the 
branching ratio n, defined as the mean aftershock 
number triggered per event, and on the sign of 9, three 
different regimes for the seismic rate N(t) are found: 



1. For n < 1 (sub-critical regime), we recover the 
results of [Sornette and Sornette, 1999], i.e. we 
find a crossover from an Omori exponent p = 
1 - 6 for t < t* to p = 1 + 6 for * > t*. 

2. For n > 1 and 9 > (super-critical regime), we 
find a transition from an Omori decay law with 
exponent p = 1 — 9 to an explosive exponential 
increase of the seismicity rate. 

3. In the case 9 < 0, we find a transition from 
an Omori law with exponent 1 — \9\ similar to 
the local law, to an exponential increase at large 
times, with a crossover time r different from the 
characteristic time t* found in the case 9 > 0. 

As we show below, these results can rational- 
ize many properties of aftershock and foreshock se- 
quences. 

The model 

We assume that a given event (the "mother") of 
magnitude rrii > too occurring at time ti gives birth 
to other events ( "daughters" ) in the time interval be- 
tween t and t + dt at the rate 

0m, (t - U) = ^_ t+c y +e H(t - U) H{ mi - m ), 

(1) 

where H is the Heaviside function: H(t — ti) = for 
t < U and 1 otherwise, mo is a lower bound magnitude 
below which no daughter is triggered. 

This temporal power law decay is Omori's law for 
the rate of aftershocks following a main shock, albeit 
with the modification that we do not specify that af- 
tershocks (daughter earthquakes) have to be smaller 
than the triggering event (mother earthquake). The 
exponential term io Q ( m_m °) describes the fact that 
the larger the magnitude m of the mother event, the 
larger is the number of daughters. The exponent 
p = 1 + 9 of the "local" Omori's law has no reason a 
priori to be the same as the one measured macroscop- 
ically which is usually found between 0.8 and 1.2 with 
an often quoted median value 1. This is in fact the 
question we address: assuming the form (1) for the 
"local" Omori's law, is the global Omori's law still a 
power law and, if yes, how does its exponent depend 
on pi What are the possible regimes of aftershocks 
as a function of the parameters of the model? 

This model can be extended to describe the spatio- 
temporal distribution of seismic activity. Following 
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Kagan and Knopoff [1981], we can introduce a spatial 
dependence in (1) of the form 



4>mi(t-ti,r-ri) 



(t-u + cy+o 



p{r-fi) H{t-U) H(m t 
(2) 



where p(r — fl) describes the probability distribution 
for an earthquake occurring at position fi to trigger 
an event an position f. This term takes into account 
the spatial dependence of the stress induced by an 
earthquake, and enable us to model the spatial dis- 
tribution of aftershocks clustered close to the main- 
shock. In this paper, we restrict our analysis to the 
temporal ETAS model without spatial dependence 
because we are mainly interested in describing the 
temporal evolution of seismic activity. The complete 
model with both spatial and temporal dependence (2) 
has been studied in [Helmstetter and Sornette, 2002] 
to derive the joint probability distribution of the times 
and locations of aftershocks including the whole cas- 
cade of secondary aftershocks. When integrating the 
rate of aftershocks calculated for the spatio-temporal 
ETAS model over the whole space, we recover the 
results given in this paper for the temporal ETAS 
model. Therefore, the results given here for the tem- 
poral ETAS model can be compared with real after- 
shock sequences when using all aftershocks whatever 
their distance from the mainshock. 

The model (1) is a branching process because each 
daughter has only one mother and not several, as 
shown in Figure 1. As we said in the introduc- 
tion, this "mean-field" assumption simplifies consid- 
erably the complexity of the process and allows for 
an analytical solution that we shall derive in the se- 
quel. The key parameter is the average number n of 
daughter-earthquakes created per mother-event. As- 
suming that the distribution P{m) of earthquake sizes 
expressed in magnitudes m follows the Gutcnberg- 
Richter distribution P(m) = b ln(10) \Q-H™-™o) ^ 
the integral of <p m {t) over time and over all magni- 
tudes m> m gives 



r+oo r + oc 
= dt 

JO J mo 



dmP(m) <j) m (t) = n a 



dt 



where 



n = 



K 



o (t + l) 1+0 ' 
(3) 

(4) 



c" b — a 

which is finite for b > a. Three cases are analyzed 
below: n < l,n = 1 and n > 1. The case n — 1 cor- 
responds to an average conservation of the number 
of events and can be associated with a brittle elas- 



tic crust without dissipation. The "dissipative" case 
n < 1 can be interpreted as corresponding to a crust 
possessing a visco-elastic component and/or a partial 
_ m cpupling with a lower ductile layer, such that a part 
of 'the energy is released aseismically. The case n > 1 
corresponds to a process in which an earthquake se- 
quence triggers an in-flow of energy from surrounding 
regions that may lead to a local sclf-cxciting ampli- 
fication. It can also correspond to a coupling with 
other non-mechanical modes of energy storage, such 
as proposed in [Sornette, 2000b; Viljocn ct al., 2002] 
which can be triggered by an event and feed the en- 
suing earthquake sequence for a while. Of course, the 
super-critical process can only be transient and has 
to cross-over to another regime. 

The case b < a requires a special attention. In ab- 
sence of truncation or cut-off, it leads to a finite-time 
singularity due to the interplay between long-memory 
and extreme fluctuations [Sornette and Helmstetter, 
2001]. However, it is more common to introduce a 
truncation or roll-off of the Gutenberg-Richter law 
at an upper magnitude. We can for example use a 
Gamma distribution of energies, which is a power- 
law distribution tapered by an exponential tail. In 
this case, the branching ratio has been calculated by 
Kagan [1991] and is given by the approximate an- 
alytical expression valid for a corner magnitude m c 
significantly larger than m , 

_ K b l0Hmc-mo) _ 10«K-">o) 
"° ~ 1? ~b~~C^ 10Kmc-mo) _ 1 ' ^) 

For a corner magnitude m c 3> mo, and for a < b, 
we recover the expression (4) for no obtained for the 
Gutenberg-Richter distribution without roll-off. 

Note that n is defined as the average over all main- 
shock magnitudes of the mean number of events trig- 
gered by a mainshock. It is thus grossly misleading 
to think of the branching ratio as giving the num- 
ber of daughters to a given earthquake, because this 
number is extremely sensitive to the specific value of 
its magnitude. Indeed, the number of aftershocks to 
a given mainshock increases exponentially with the 
mainshock magnitude as given by (1), so that large 
earthquakes will have many more aftershocks than 
small earthquakes. ^From (1) and (3), we can calcu- 
late the mean number of aftershocks N(M) triggered 
directly by a mainshock of magnitude M 

(b-a). 



N(M) = n- 



(6) 



As an example, take a = 0.8, b = 1, m = and 
n = 1. Then, a mainshock of magnitude M = 7 will 
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have on average 80000 direct aftershocks, compared 
to only 2000 direct aftershocks for an earthquake of 
magnitude M = 5 and less than 0.2 aftershocks for 
an earthquake of magnitude M = 0. 

When 9 > 0, J °° ( t+ fy+e = 1/9 and the branch- 
ing ratio n = n /9 is finite. In this regime, n is an 
increasing function of the rate K and a decreasing 
function of 9, c and b — a. 

Even for b > a and 9 > 0, the average number of 
daughters per mother can be larger than one: n > 1. 
This regime corresponds to the super-critical regime 
of branching processes [Harris, 1963; Sornette, 2000a] 
in which the total number of events grows on average 
exponentially with time. If n < 1, there is less than 
one earthquake triggered per earthquake on average. 
This is the sub-critical regime in which the number 
of events following the first main shock decays even- 
tually to zero. The critical case n = 1 is at the bor- 
derline between the two regimes. In this case, there 
is exactly one earthquake on average triggered per 
earthquake and the process is exactly at the critical 
point between death on the long run and exponential 
proliferation. 

There is another scenario, occurring for 9 < 0, 
in which the seismicity blows up exponentially with 
time. In this case, the integral J °° ( t+ f*i + o becomes 
unbounded. In principle, n becomes infinite: this 
does not invalidate the ETAS model per se. It only re- 
flects the fact that the calculation of an average num- 
ber of daughters per mother has become meaningless 
because of the anomalously slow decay of the kernel 
(j>(t) . This mechanism is reminiscent of that leading to 
anomalous diffusion and to aging in quenched random 
systems and spinglasses (see [Sorncttc, 2000a] for an 
introduction) . As in these systems, any estimation of 
the averages depend on the time scale of study: due 
to the extremely slow decay of (f>(t), the number of 
daughters created beyond any time t far exceeds the 
number of daughters created up to time t. Notwith- 
standing the decay, its cumulative effect creates this 
dominance of the far future. This regime is the op- 
posite of the situation where 9 > where most of the 
daughters are created at relatively early times. Since 
the number of daughters born up to time t is an un- 
bounded increasing function of t, it is intuitively ap- 
pealing, as we show in the appendix, that this regime 
should be similar to the super-critical regime n > 1 
discussed above in the case 9 > 0. 

Until now, we have discussed three issues related 
to the convergence of the ETAS sequences: (i) the 



condition 9 > ensures convergence at large times; 
(ii) the convergence at short times is obtained by the 
introduction of the regularization constant c in the 
generalized Omori's law; (iii) the condition a < b is 
a necessary condition for the finitcness of the number 
of daughters. Finally, we should stress the role of the 
"ultra-violet" cut-off mo on the magnitudes. In the 
ETAS model, only earthquakes of magnitude m > mo 
are allowed to give birth to aftershocks, while events 
of smaller magnitudes are lost for the epidemic dy- 
namics. If such a cut-off is not introduced and no 
cut-off is put on the Gutenberg-Richter toward small 
magnitudes, the dynamics becomes completely domi- 
nated by the swarms of very tiny earthquakes, which 
individually has very low probability to generate af- 
tershocks but become so numerous that their collec- 
tive effect becomes overwhelming in the dynamics. 
We would thus have the unphysical situation in which 
a magnitude 7 or 8 earthquake may be triggered by 
tiny earthquakes of magnitudes —2 or less. We stress 
that the introduction of such a cut-off mo is a simple 
way to prevent such a situation to occur, but it does 
not mean that small earthquakes of magnitude below 
mo do not have their own aftershocks. It only means 
that such small earthquakes create aftershocks that 
can not participate in the epidemic process leading 
to significantly larger earthquakes; these small earth- 
quakes live their separate life. This is why they are 
not registered by the ETAS model. This formula- 
tion is of course only an end-member of many possi- 
ble regularization procedures, which are well-known 
to be an ubiquitous requisite in mechanical models of 
rupture. An improvement of the ETAS model would 
be for instance to replace this abrupt cut-off mo by 
introducing a roll-off in the Gutenberg-Richter law 
for the aftershocks with a characteristic corner mag- 
nitude decreasing with the magnitude of the mother 
earthquake. This and other schemes will not be ex- 
plored here, as we want to analyze the simplest ver- 
sion possible. 

We now describe briefly the connection with pre- 
vious works in the mathematical statistics literature. 
As we said above, the model (1) belongs to the gen- 
eral class of branching models [Moyal, 1962; Har- 
ris, 1963]. The elementary results on the stability 
of the process, such as the condition n < 1, have 
been known for many years, and go back to the ori- 
gin of the ETAS model as a special case (for dis- 
crete magnitudes) or extension (for continuous mag- 
nitudes) of the class of "mutually exciting point pro- 
cesses" introduced in [Hawkes, 1971; 1972; Hawkes 
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and Adamapoulos, 1973]. A convenient mathemati- 
cal overview is in Chapter 5 of Vere- Jones and Da- 
ley [1988], especially Example 5.5(a) and associated 
exercises 5.5.2-5.5.6. For the ETAS model, the equa- 
tions governing the probability generating functional, 
the probability of extinction within a given number 
of generations, the expectation measure for the total 
population, the second factorial moment (related to 
the covariance of the population) and their Fourier 
transform can be derived as special cases of results 
summarized there. In particular, the process initiated 
with a single event at the origin corresponds to the 
total progeny process for a general branching process 
model with time-magnitude state space and a single 
ancestor at time t = 0; Exercise 5.5.6 gives the equa- 
tions of the above cited variables for the case of fixed 
magnitudes (i.e., a = 0). This direct probabilistic 
analysis in terms of generating functions effectively 
replaces the Wiener- Hopf theory in the present paper 
and mentioned also in [Hawkes, 1971; 1972; Hawkes 
and Oakes, 1974]. However, there is not explicit solu- 
tions given to these equations and there is no discus- 
sion of the change of regime from an effective Omori's 
law \jt x ~ 6 at early times to l/t 1+e at long times, nor 
mention of the interesting super-critical done 
in the present work. 

Hawkes [1971; 1972] and Hawkes and Adamapou- 
los [1973] use what is in effect an ETAS model with 
an exponential "bare" Omori's law rather than the 
power law l/(i + c) 1+e defined in (1). Hawkes and 
Adamapoulos [1973] use it in an early study of earth- 
quake data. The introduction of magnitudes is simi- 
lar to the introduction of a marked process associated 
with a single point process [Hawkes, 1972]; however, 
the impact of magnitudes on the seismicity rate is as- 
sumed to be linear in [Hawkes, 1972] while it is mul- 
tiplicative in the ETAS model. Our derivation pre- 
sented in the appendix of the solution of the ETAS 
model for the mean rate of earthquakes in terms of 
its Laplace transform recovers previous results. For 
instance, equation (17) in [Hawkes and Oakes, 1974] 
is the same as our equation (A6) in our Appendix 
(up to a factor (3 stemming from taking the cumu- 
lative number in [Hawkes and Oakes, 1974]). The 
key factor Q{f3) in (A7) corresponds to the quantity 
Gi(0) in equation (5) of [Hawkes, 1972]. The link 
between Hawkes' "mutually exciting point processes" 
and branching processes was made explicit in [Hawkes 
and Oakes, 1974]. 

Some average properties of the ETAS model have 
been derived in the Master thesis of P.A. Ramsc- 



laar (1990). Specifically, using the theory of Markov 
processes applied to branching processes, Ramselaar 
(1990) proves that, in the supercritical regime n > 1 
(where n is the average branching ratio defined in (5)), 
the average number of aftershocks stemming from 
a common ancestor grows exponentially as ~ e*'* 
where t* is the solution of nR(c/t*) = 1 and the 
function R is defined in (A9). The solution of this 
equation nR(c/t*) = 1 for t* is the same as our t* 
given by (12) and the exponential growth of Ramse- 
laar is therefore the same as our result (17). We add 
on this asymptotic result, which is valid only at large 
times, by exhibiting the solution for the aftershock 
decay at early times. In addition, contrary to the in- 
correct claim of Ramselaar (1990) that "the Ogata 
earthquake process is critical or supercritical but is 
never subcritical," we demonstrate that the subcriti- 
cal regime exhibits a rich phenomenology 

Analytical solution 

We analyze the case where there is an origin of 
time t = at which we start recording the rate of 
earthquakes, assuming that the largest earthquake of 
all has just occurred at t = and somehow reset 
the clock. In the following calculation, we will forget 
about the effect of events preceding the one at t = 
and count aftershocks that are created only by this 
main shock. 

Let us call N m (t) the rate of seismicity at time 
t and at magnitude m, that is, N m (t)dtdm is the 
number of events in the time/magnitude interval 
dt x dm. We define its expectation X m (t)dtdm = 
E[N m (t) dtdm], as the mean number of earthquakes 
occurring between t and t + dt of magnitude between 
to and to + dm. X m (t) is the solution of a self- 
consistency equation that formalizes mathematically 
the following process : an earthquake may trigger af- 
tershocks; these aftershocks may trigger their own af- 
tershocks, and so on. The rate of seismicity at a given 
time t is the result of this cascade process. The self- 
consistency equation that sums up this cascade reads 

ft 



X m (t) = E[N m (t)] 



E 



- f 

J m 



poo p 

/ dm' / 

J mo J - 

7. 



dr (j) m ,(t-T)P{m)N n 



dm 



dm' 



dT <j) m ,{t-T)E[N m ,{T)] 



dr <j> m ,(t - r)A m /(r) . 



If there is an external source S(t,m), it should be 
added to the right-hand-side of (9). 
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The mean instantaneous rate A TO (t) at time t is 
the sum over all induced rates from all earthquakes 
of all possible magnitudes that occurred at all pre- 
vious times. The rate of events at time t induced 
per earthquake that occurred at an earlier time r 
with magnitude m! is equal to (p m i (t — r) . The term 
P(m) is the probability that an event triggered by 
an earthquake of magnitude m' is of magnitude m. 
We assume that this probability is independent of 
the magnitude of the mother-earthquake and is noth- 
ing but the Gutcnbcrg-Richter law. This hypothesis 
can be easily relaxed if needed and P(m) can be gen- 
eralized into P(m\m') giving the probability that a 
daughter-earthquake is of magnitude m conditioned 
on the value vn! of the magnitude of the mother- 
earthquake. However, we do not pursue here this pos- 
sibility as this hypothesis seems well-founded empiri- 
cally [Ranalli, 1969; Knopoff et al, 1982]. The term 
S(t, m) is an external source which is determined by 
the physical process. We consider the case where a 
great earthquake occurs at the origin of time t = 
with magnitude M. In this case, the external source 
term is 

S(t, m) — S(t) S(m - M) , (10) 

where 8 is the Dirac distribution. Other arbitrary 
source functions can be chosen. 

By construction of the kernel (1), the integrals over 
the magnitudes and over time factorize, which implies 
that the solution for X m (t) can be searched as 



X m (t) = P(m)X(t) , 



(11) 



which makes explicit the separation of the variables 
magnitude and time. We stress that (11) is not an 
assumption: it is the structure of the solution based 
on the assumption made at the level of equations (1) 
and (9). Therefore, the ETAS model assumes that the 
time response and the magnitude response are inde- 
pendent. In reality and more generality, wc can envi- 
sion that the rate of activation of new earthquakes will 
depend on 1) the magnitude of the "mother" (which 
the ETAS model takes into account multiplicatively 
in (1)), 2) on the magnitude of the daughter (which 
is neglected in the ETAS model) and 3) on the time 
since the mother was born. Rather than having a 
very general kernel combining these three parameters 
nonlinearly, equations (1) and (9) are based on an 
hypothesis of independence between these different 
factors that allows us to factorize them, leading to 
(11). 

The problem is then to determine the functional 
form of A(t), assuming that <j) is given by (1). The in- 



tegral equation (9) is a Wiener-Hopf integral equation 
[Feller, 1971]. It is well-known [Feller, 1971; Morse 
and Feshbach, 1953] that, if <P(t) decays no slower 
than an exponential, then X(t) has an exponential 
tail X(t) ~ exp[— rt] for large t with r solution of 
J 4>( x ) exp[rx] dx = 1. This result implies that a 
global Omori's law cannot be obtained by the epi- 
demic ETAS branching model with, for instance, lo- 
cal exponential relaxation rates. In the present case, 
</>(r) decays much slower than an exponential and a 
different analysis is called for that we now present. 
The solution of (9) is derived in the Appendix and is 
summarized in the following sections. For the sequel, 
it is useful to define the characteristic time 



t =ct 



0) 



l-n 



(12) 



where T(x) is the Gamma function: T(z) = duu z 1 e 
which is nothing but (z — 1)! for positive integers z. 

The sub-critical regime n < 1 and 9 > 

An approximation is made in the analytical solu- 
tion so that the results presented below are only valid 
for t > c. 

We define the parameter So that describes the ex- 
ternal source term 



5 = ^ ^) 10 a(M-mo) _ 



Two cases must be distinguished. 
• For c<f«f, we get 



(13) 



A t <*.(t) 



So 



t* 



r(0)(l-n) t 1 - 8 

For t > t* , we obtain 

S t* e 



X t >t*(t) 



r(0)(l-n) t 1+e 



for c < t < t* . 

(14) 

fort>r . (15) 



We verify the self-consistency of the two solutions 
X t >t*(t) and X t<t *(t) by checking that \ t >r(t*) = 
Xt<t*(t*)- In other words, t* is indeed the transition 
time at which the "short-time" regime A t<t » (t) crosses 
over to the "long-time" regime X t>t *(t). 

The full expression of X(t) valid at all times t ^> c 
is given by 



Mt) = j 



So t*~ 



(-i) A 



„ t i-e Z_.v -> T((k + 1)0) 



(16) 
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Expression (16) provides the solution that describes 
the cross-over from the \jt x ~ 6 Omori's law (14) at 
early times to the l/t 1+e Omori's law (15) at large 

times. The series Z)fclo(~-O fc r*(fc+i)fl) ^ s a ser i es rep- 
resentation of the special Fox function [Glockle and 
Nonnenmacher, 1993] (see the Appendix for details). 

The ETAS model has been simulated numerically 
using the algorithm described in [Ogata, 1998, 1999]. 
Starting with a large event of magnitude M at time 
t = 0, events are then simulated sequentially. After 
each event, we calculate the conditional intensity X{t) 
defined by 

_ ^ K io Q ( m *-™o) 

<<)= 6*,c-«+<> I+ * 

where t is the time of the last event and t% and rm 
are the times and magnitudes of all preceding events 
that occurred at time U < t. The time of the follow- 
ing event is then determined according to the non- 
stationary Poisson process of conditional intensity 
X(t), and its magnitude is chosen in a Gutenberg- 
Richter distribution with parameter b. Theses sim- 
ulations are compared to the theoretical predictions 
in Figure 2, which shows the aftershock seismic rate 
A(i) in the sub-critical regime triggered by a main 
event of M = 6.8, for the parameters K = 0.024 
(constant in (1)), the threshold mo = for aftershock 
triggering, c = 0.001, a — 0.5, a &-valuc b = 1.0 and 
9 = 0.2 (corresponding to a local Omori's exponent 
p = 1.2). These parameters lead to a branching ratio 
n = 0.95 (equation (3)) and a characteristic cross-over 
time t* = 4500 (equation (12)). The noisy black line 
represents the seismicity rate obtained for the syn- 
thetic catalog. The local Omori law with exponent 
p = 1 + 6 = 1. 2 is shown for reference as the dot- 
ted line. The analytical solution (16) is shown as the 
thick line. The two dashed lines represent the approx- 
imation solutions (14) for t < t* and (15) for t > t*. 

The super-critical regime n > 1 and 9 > 

From the definition of the branching ETAS model 
for n > 1, it is clear that the number of events X(t) 
blows up exponentially for large times as n — 1 to a 
power proportional to the number t of generations. 
We shall show below that the rate of the exponen- 
tial growth can be calculated explicitly, which yields 
A(t) ~ e*/**, where t* has been defined in (12). How- 
ever, there is an interesting early and intermediate 
time regime in the situation where a great event of 
magnitude M has just occurred at t = 0. In this 



case, the total seismicity is the result of two com- 
peting effects: (1) the total seismicity tends to de- 
cay according to the Omori's law governing the rate 
of daughter-earthquakes triggered by the great event; 
(2) since each daughter may in turn trigger grand- 
daughters, grand-daughters may trigger grand-grand 
daughters and so on with a number n > 1 of chil- 
dren per parent, the induced seismicity will eventually 
blow up exponentially. However, before blowing up, 
one can expect that seismicity will first decay because 
it is mainly controlled by the large rate ~ io a ( M ~ m °) 
directly induced by the great earthquake which de- 
cays according to its "local" Omori's law. This decay 
will be progressively perturbed by the proliferation of 
daughters of daughters of ... and will cross-over to 
the explosive exponential regime. 

At early times cCf<f, the early decay rate of 
aftershocks is the same « (S /T(6)(n-1)) {t*~ e /t 1 - ) 
as for the sub-critical regime (14) (see the Appendix). 
However, as time increases, the Appendix shows that 
the decay of aftershock activity can be represented 
as a power law with an effective apparent expo- 
nent # app > 9 increasing progressively with time. 
The seismic rate will thus decay approximately as 
~ l/f 1_e *pp(*). Quantitatively, the large time behav- 
ior is (see the Appendix) 

X(t) ~ S \ - e*'*' (17) 
w (n-l)t*6 K ' 

exhibiting an exponential growth at large times. Ex- 
pression (12) shows that 1/t* ~ |1 — n\t . Thus, as ex- 
pected, the exponential growth disappears as n — > 1 + . 
The full expression of X(t) valid at times t ^> c is 

m -Jn-l)¥^^Y{(k + l)9) (18) 

Expression (18) provides the solution that describes 
the cross-over from the l/t 1 ^ 9 Omori's law at early 
times (14) to the exponential growth (17) at large 
times. 

Figure 3 tests these predictions by comparing them 
with direct numerical simulation of the ETAS model, 
in the case of a main shock of magnitude M = 6. The 
parameters of the synthetic catalog are K = 0.024 
(constant in (1)), the threshold mo = for aftershock 
triggering, c = 0.001, a — 0.5, a b- value b = 0.75 and 
9 = 0.2 (corresponding to a local Omori's exponent 
p = 1.2). These parameters lead to a branching ratio 
n = 1.43 (equation (3)) and a characteristic cross- 
over time t* = 0.85 (equation (12)). The noisy black 
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line represents the seismicity rate obtained for the 
synthetic catalog. The local Omori law with exponent 
p = 1 + = 1. 2 is shown for reference as the dotted 
line. The analytical solution (18) is shown as the 
thick line. The two dashed lines correspond to the 
approximative analytical solutions (14) and (17). At 
early times c < t < t* , the decay of N(t) is initially 
close to the prediction (14). For t > t* , we observe 
that the analytical equation (18) is very close to the 
exponential solution (17). 

Case 9 < corresponding to a local Omori's 
law exponent p < 1 

We have already remarked that, in this case, the in- 
tegral Jq 00 (j. + fy +e in the definition (3) of the branch- 
ing ratio n becomes unbounded: the number of daugh- 
ters created beyond any time t far exceeds the number 
of daughters created up to time t. 

The appendix shows that the general equation (9) 
still holds and the general derivation starting with 
(Al) up to (A6) still applies. 

Similarly to the super-critical case n > 1 of the 
regime 9 > 0, we find a crossover from a power-law 
decay at early times to an exponential increase of the 
seismicity rate at large times. The characteristic time 
t that marks the transition between these two regimes 
is given by 



( n T(\e\) \ 

\ 1 + wJ 



(19) 



In contrast with the case 8 > 0, the early time 
behavior (i.e., c<Kt) of the global decay law in 
the case 9 < is similar to the local Omori law : 



A(t) 



(i + ffl)r(|0|) 



(20) 



Similarly to the super-critical case n > 1 of the 
regime 9 > 0, the long time dependence of the regime 
9 < is controlled by a simple pole 1/r leading to a 
long-time seismicity growing exponentially 



X(t) 



So 



(21) 



This result is in agreement with the fact that the num- 
ber of daughters born up to time t is an unbounded 
increasing function of t, and we should thus recover a 
regime similar to the super-critical case of 9 > 0. 



The full expression of X(t) valid at times t > c is 
Mt) So 1 V {t/T)m (22) 

Expression (22) provides the solution that describes 
the cross-over from the local Omori law 1/i 1 "' 6 at 
early times to the exponential growth at large times. 

Figure 4 compares these predictions to a direct nu- 
merical simulation of the ETAS model, in the case of a 
main shock of magnitude M = 7. The parameters of 
the synthetic catalog are K = 0.02, m = 0, c = 0.01, 
a = 0.5, 6=1 and 9 = —0.1 (corresponding to a local 
Omori's exponent p = 0.9). These parameters lead 
to a characteristic cross-over time r = 10 5 (equation 

(19) ). The noisy black line represents the seismic- 
ity rate obtained for the synthetic catalog. The local 
Omori law with exponent p = 1 + 9 = 0.9 is shown for 
reference as the dotted line. The analytical solution 
(22) is shown as the thick line. The two dashed lines 
correspond to the approximative analytical solutions 

(20) and (21). At early times c < t < t, the decay of 
A(i) is initially close to the prediction (20). For t > r, 
we observe that the analytical equation (22) is very 
close to the exponential solution (21). 

Discussion 

Assuming that each event triggers aftershock se- 
quences according to the local Omori law with expo- 
nent 1 + 6, we have shown that the decay law of the 
global aftershock sequence is different from the local 
one. Depending on the branching ratio n, which is a 
function of all ETAS parameters, we find two differ- 
ent regimes, the sub-critical regime for n < 1 and the 
super-critical regime for n > 1 and 9 > 0. For the 
two regimes in the case 9 > 0, a characteristic time 
t* , function of c, n and 9, appears in the global decay 
law \(t) and marks the transition between the early 
time behavior and the large time behavior. In the 
sub-critical regime (n < 1), the global decay law is 
composed of two power laws. At early times (t < t*), 
A(t) decays like t~ 1+e . At large times (t > t*) the 
global decay law recovers the local law N(t) ~ t^ 1 ^ 6 . 
In the super-critical regime (n > 1 and 9 > 0), the 
early times decay law is similar to that of the sub- 
critical regime, and the seismicity rate increases ex- 
ponentially for large times. The case 9 < leads to 
an infinite n-value, due to the slow decay with time 
of the local Omori law. In this case, we find a transi- 
tion from an Omori law with exponent 1 — \9\ similar 
to the local law, to an exponential increase at large 
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times, with a crossover time r different from the char- 
acteristic time t* found in the case 9 > 0. Thus, the 
Omori law is only an approximation of the global de- 
cay law valid for some time periods and parameter 
values. The value of the local Omori exponent p = 1 
is the only one for which the local and the global de- 
cay rate are similar, and are both power-laws without 
any characteristic time. For small n, t* is very small 
so that in real data we should observe only the be- 
havior t > t* characteristic of large times. The global 
decay law then appears similar to the local Omori 
law. On the contrary, for n close to 1, t* is very large 
by comparison with the time period available in real 
data, and we should observe only the power-law be- 
havior A(i) ~ t~ 1+e characteristic of early times, with 
a global p- value smaller than the local one. Changing 
n thus provides an important source of variability of 
the exponent p. 

Estimation of n and t* in earthquake data 

In real earthquake data, it is possible to evaluate 
the branching value n in order to determine if the seis- 
mic activity is either in the sub- or the super-critical 
regime. The values of n and t* can be evaluated from 
equations (3) and (12) as a function of the ETAS pa- 
rameters b, p = 1 + 9, c, K and a. The parameters 
of the ETAS model and their standard error can be 
inverted from seismicity data (time and magnitudes 
of each event) using a maximum likelihood method 
[Ogata, 1988]. We now discuss the range of the dif- 
ferent parameters obtained from such inversion pro- 
cedure. 

• The parameter a is found to vary between 0.35 
to 1.7, and is often close to 0.5 [Ogata, 1989, 
1992; Quo and Ogata, 1997]. An a- value of 0.5 
means that a mainshock of magnitude M will 
have on average 10 times more aftershocks than 
a mainshock of magnitude M — 2, independently 
of M. Note that our definition of a is slightly 
different from that used by Ogata and we have 
divided his a-values by ln(10) to compare with 
our definition. 

For some seismicity sequences, Ogata [1989, 
1992] and Quo and Ogata [1997] found a > b. 
According to (3), this leads to an infinite re- 
value if we use a Gutenberg-Richter magni- 
tude distribution. As we said, a truncation of 
the magnitude distribution is needed to obtain 
a physically meaningful finite n-value because 
the seismicity rate is controlled by the largest 



events. 

A large a-value can be associated with seismic 
activity called "swarms" , while a small a-value 
is observed for aftershock sequences with a sin- 
gle mainshock and no significant secondary af- 
tershock sequences [Ogata, 1992, 2001]. 

• The parameter c is usually found to be of the or- 
der of one hour [Utsu et al, 1995]. In practice, 
the evaluation of c is hindered by the incom- 
pleteness of earthquake catalogs just after the 
occurrence of the mainshock, due to overlap- 
ping aftershocks on the seismograms. A large c 
is often an artifact of a change of the detection 
threshold. Notwithstanding these limitations, 
well-determined non-zero c-value have been ob- 
tained for some aftershocks sequences [Utsu et 
al, 1995]. Note that a non-zero c is required for 
the aftershocks rate to be finite just at the time 
of the mainshock. 

• The "local" p-value, equal to 1 + 9, describes 
the decay law of the aftershock sequence trig- 
gered by a single earthquake. The local Omori 
law is the law <p(t) obtained by inverting the 
ETAS model on the data. The "global" p- 
value describes the decay law of the whole af- 
tershock sequence, composed of all secondary 
aftershocks triggered by each aftershock. The 
global Omori law is the law </>(i) fitted directly 
on the data, without taking into account the 
hierarchical structure of branching of the ETAS 
model. We have shown that the Omori law is 
only an approximation of the global decay law, 
so that in the subcritical regime the global p- 
value will change from 1 — 9 at early times to 
1+9 at large times. [Guo and Ogata, 1997] mea- 
sured both the local and global p- values for 34 
aftershocks sequences in Japan, and found that 
the local p- value is usually slightly larger than 
the global p- value [Guo and Ogata, 1997]. This 
is in agreement with our prediction when iden- 
tifying the local p- value with 1 + 9 (recovered at 
large times) and the global p-value with 1 — 9 
found at early times. Guo and Ogata [1997] 
and Ogata [1992, 1998, 2001] found a local p- 
value smaller than one for some aftershocks se- 
quences in Japan. Within the confine of the 
ETAS model, this corresponds to the case 9 < 
discussed above and in the appendix. 

• The parameter K measures the rate of after- 
shocks triggered by each earthquake, indepen- 
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dcntly of its magnitude. Recall that the branch- 
ing ratio n is proportional to K. It is usually 
found of the order of K w 0.02 [Ogata, 1989, 
1992; Guo and Ogata, 1997], but large varia- 
tions of -fT-value from 0.001 to 5 are reported 
by Ogata [1992]. 

• The parameter /i measures the background seis- 
micity rate that is supposed to arise from the 
tectonic loading, /i ~ for an aftershock se- 
quence triggered by a single mainshock. This 
parameter as no influence on the branching ra- 
tio n. In real catalogs, the background scismic- 
ity only accounts for a small part of the seismic 
activity. 

We have computed the branching ratio n and the 
cross-over time t* from the ETAS parameters mea- 
sured by Ogata [1989, 1992] for several seismicity se- 
quences in Japan and elsewhere. The ETAS param- 
eters and the n and t* values are given in Table 1. 
When the &-value is not given in the text, we have 
computed n and t* assuming a b- value is equal to 1. 
We find that the n-value is either smaller or larger 
than 1. This means that the seismicity can be inter- 
preted to be either in the sub- or in the super-critical 
regime. An infinite n-value is found if the local p- 
value is smaller than one (8 < 0) or if the a-value 
is larger than the &-value. For the same area, the 
ETAS parameters and the n and t* values are found 
to vary in time, sometimes changing from the sub- to 
the super-critical regime. The characteristic time t* 
shows large spatial and temporal variability, ranging 
from 0.4 days to 10 22 days. Large t* values are re- 
lated to a branching ratio n close to one, i.e., close 
to the critical point n = 1. The ETAS model thus 
provides a picture of seismicity in which sub-critical 
and super-critical regimes are alternating in an inter- 
mittent fashion. As we shall argue, the determination 
of the regime may provide important clues and quan- 
titative tools for prediction. 

Implications of the ETAS model in the 
sub-critical regime n < 1 

In the sub-critical regime, the ETAS model can ex- 
plain many of the departures of the global aftershock 
decay law from a pure Omori law. 

The ETAS model contains by definition (and thus 
"explains") the secondary aftershock sequences trig- 
gered by the largest aftershocks that are often ob- 
served [Correig et at, 1997; Guo and Ogata, 1997; 
Simeonova and Solakov, 1999; Ogata, 2001]. In the 



ETAS model, the fact that secondary aftershock se- 
quences of large aftershocks can stand out above the 
overall background aftershock seismicity results from 
the factor io«( m «- m ») in (1). 

Our analytical results may rationalize why some al- 
ternative models of aftershock decay work better than 
the simple modified Omori law. In the sub-critical 
regime, we predict an increase of the apparent global 
p- value from 1 — 8 at early times to 1+6* at large times. 
To our knowledge, this change of exponent has never 
been observed. This change of power law may be 
approximated by the stretched exponential function 
proposed by [Kissinger, 1993; Gross and Kisslinger, 
1994] to fit aftershocks sequences. In the stretched 
exponential model, the rate of aftershocks X(t) is de- 
fined by 

\(t) = K ^9- 1 e -( t / t '') , , (23) 

where q, K and to are constants. At early times, this 
function decays as a power law l/t l ~ q with appar- 
ent Omori's exponent 1 — q. For times larger than 
the relaxation time to, the seismicity rate decays ex- 
ponentially in the argument (t/to) q . For q < 1, this 
decay is much slower than exponential and can be 
accounted for by an apparent power law with larger 
exponent. Figure 5 compares the stretched expo- 
nential function with the analytical solution of the 
ETAS model (16) with parameters t* = t and 8 = q, 
and with the Omori law of exponent p = 1 — q. 
These three laws have the same power-law behav- 
ior at early times, and then both the stretched expo- 
nential and the analytical solution (16) decay faster 
than the Omori law at large times. The fact that 
it is very difficult to distinguish the decay laws de- 
scribed by power laws and by stretched exponential 
has been illustrated in [Laherrere and Sornette, 1998] 
in many examples including earthquake size and fault 
length distributions. Kissinger [1993] and Gross and 
Kisslinger, 1994] compared this function to the mod- 
ified Omori law X(t) = K (t + c)~ p for several after- 
shocks sequences in southern California. They found 
that the stretched exponential fit often works better 
for the sequences with a small p-valuc or a large q- 
value, indicative of a slow decay for small times. This 
is in agreement with our result that in the sub-critical 
regime a slowly decaying aftershock sequence (global 
p-value smaller than one) will then cross-over to a 
more rapid decay for time larger than t* . The relax- 
ation time to ranges between 2 days and 380 days for 
the sequences that are better fitted by the stretched 
exponential [Kissinger, 1993]. This parameter is anal- 
ogous to t* found in our model, because these two 
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parameters define the transition from the early time 
power-law decay to another faster decaying behavior 
for large times. To further validate our results, these 
aftershocks sequences should be fitted using equation 
(16) to compare our results with the stretched expo- 
nential function and determine if the transformation 
of the early time power law decay is better fitted by a 
stretched exponential fall-off or an increase in the ap- 
parent Omori exponent from 1 — 9 to 1 + 9 as predicted 
by our results. 

The ETAS model can also rationalize some cor- 
relations found empirically between seismicity pa- 
rameters. It may explain the rather large variabil- 
ity of the global empirical p-valuc. Guo and Ogata 
[1995] have reported a positive correlation between 
the Gutenberg- Richter b- value and the p- value (expo- 
nent of the global Omori law) for several aftershocks 
sequences in Japan. A similar correlation has also 
been found by [Kisslinger and Jones, 1991] for sev- 
eral aftershock sequences in southern California, but 
this correlation was detectable only if the earthquake 
sequences were separated into thrust and strike slip 
events. This positive correlation between b and global 
p values is expected from our analysis. From equation 
(3), we see that a small b- value is associated with a 
large n value. For n ~ 1, the characteristic time t* is 
very large, so that the global aftershock rate decays as 
a power law with exponent 1 — 9 over a large time in- 
terval. For n > 1 and 9 > 0, we see an apparent global 
p- value smaller than 1 — 9 which decreases with time. 
In contrast, for large ^-values, the branching ratio n 
is small and the characteristic time t* is very small. 
In this case, only the large time behavior is observed 
with a larger exponent 1 + 9. Consequently, in the 
subcritical regime, our results predict a change of the 
global p- value from 1 — 9 for small b- value and times 
t <C t* to 1 + 9 for large b- values. There is also a pos- 
itive correlation between p-value and &-value in the 
super-critical regime. For n > 1 or 9 < 0, the global 
aftershock sequence is characterized by an apparent 
exponent p smaller than 1 — |0| which decreases with 
time. Then, we expect the apparent exponent p to 
be all the smaller, the smaller is the 6-value, because 
the characteristic times t* for 9 > or t for 9 < 
decreases with b. The variability of the global p expo- 
nent reported by Guo and Ogata [1995] and Kisslinger 
and Jones [1991] may thus be explained by a change 
of b- value and a constant local p exponent. However, 
the results of Guo and Ogata [1997] contradict this in- 
terpretation. Guo and Ogata [1997] studied the same 
aftershocks sequences than Guo and Ogata [1995] but 



they measured the local p-value of the ETAS model. 
They still found a large variability in the local p- value, 
and a positive correlation between this local p- value 
and the b- value. 

Implications of the ETAS model in the 
super-critical regime and in the case 9 < 

In the regime where the mean number of after- 
shocks per mainshock is larger than one (i.e., n > 1), 
the mean rate of aftershocks increases exponentially 
for large times. However, because of the statisti- 
cal fluctuations, the aftershock sequence has a finite 
probability to die. This probability of extinction can 
be evaluated for the simple branching model without 
time dependence [Harris, 1963]. Therefore, a branch- 
ing ratio larger than 1 does not imply necessarily that 
the number of aftershocks will be infinite. If n is not 
too large, and if the number of aftershocks is small, 
there is a significant probability that the aftershock 
sequences will die, as observed in numerical simula- 
tions of the ETAS model. If the characteristic time t* 
is very large, the aftershock sequence may not remain 
supercritical long enough for the exponential increase 
to be observed. Even if the large times exponential ac- 
celeration is rarely observed in real seismicity, it may 
explain the acceleration of the deformation before ma- 
terial failure. The early times behavior of the seis- 
mic activity preceding the exponential increase has 
also important possible implications for earthquake 
prediction, and can rationalize some empirically pro- 
posed seismic precursors, such as the low p- value [Liu, 
1984; Bowman, 1997], or the relative seismic quies- 
cence preceding large aftershocks [Matsu'ura, 1986; 
Drakatos, 2000]. 

It is widely accepted that about a third to a half of 
strong earthquakes are preceded by foreshocks [e.g., 
Jones and Molnar, 1979; Bowman and Kisslinger, 
1984; Reasenberg, 1985, 1999; Reasenberg and Jones, 
1989; Abercombie and Mori, 1996], i.e., are preceded 
by an unusual high seismicity rate for time peri- 
ods of the order of days to years, and distance up 
to hundreds kilometers. However, there is no reli- 
able method for distinguishing foreshocks from after- 
shocks. Indeed, the ETAS model makes no arbitrary 
distinctions between foreshocks, mainshocks and af- 
tershocks and describes all earthquakes with the same 
laws. While this seems a priori paradoxical, our anal- 
ysis of the ETAS model provides a useful tool for iden- 
tifying foreshocks, i.e., earthquakes that are likely to 
be followed by a larger event, from usual aftershocks 
that are seldom followed by a larger earthquake. The 
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characterization of foreshocks will be performed in 
statistical terms rather than on a single-event basis. 
In other words, we will not be able to say whether 
any specific event is a precursor. It is the ensemble 
statistics that may betray a foreshock structure. 

The crux of the method is that, when seismicity 
falls in the regime with a branching ratio n > 1, 
the corresponding earthquake sequences can be iden- 
tified as foreshocks. This is because the super-critical 
regime corresponds to an exponentially accelerating 
seismicity for times larger than t*: by a pure statis- 
tical effect, the larger number of earthquakes of any 
size will sample more and more the branch of the 
Gutenbcrg-Richter law toward large events. Thus by 
the sheer weight of numbers, larger and larger earth- 
quakes will occur as time increases. Of course, we are 
not implying any precise deterministic growth law, 
but statistically, the largest events should indeed grow 
significantly, the more so, the more within the super- 
critical regime, the larger the branching ratio n > 1. 
Conversely, this argument implies that, in the subcrit- 
ical regime, the triggered events are usual aftershocks, 
because a mainshock is unlikely to be followed by a 
larger triggered event. Foreshock sequences can thus 
be identified by evaluating the branching ratio n from 
the inversion of seismic data (times and magnitudes 
of an earthquake sequence) for the ETAS parameters. 
There is however a finite probability than a triggered 
event in the subcritical regime be larger than the trig- 
gering event, and thus the triggering event will be 
a foreshock of the triggered event. Therefore, fore- 
shocks can be observed even in the sub-critical regime, 
but they are less frequent than aftershocks. 

A note of caution is in order: the direct estima- 
tion of n and t* or t may be quite imprecise if the 
number of events is small. Based on our analysis and 
our results, the foreshock regime can be nevertheless 
identified with relatively good confidence if one as- 
sumes an upper bound for the local exponent p. Let 
us assume for instance that the local p- value is smaller 
than 1.3 (i.e., 9 < 0.3); according to our results, the 
global exponent p = 1 + 8 cannot become smaller than 
1 — 9 = 0.7 in the sub-critical regime. In contrast, in 
the supercritical regime, we have shown that the ap- 
parent exponent is smaller than or at most equal to 
1 — 9. Therefore, a measure of the global p- value yield- 
ing a value smaller than 0.7, is always associated with 
the super-critical regime. As we said above, Guo and 
Ogata [1997] and Ogata [1992, 1998, 2001] found a lo- 
cal p-value smaller than one for some aftershocks se- 
quences in Japan corresponding to the case < 0. A 



small global p- value can thus also result from a small 
local p-value. In sum, a small global p-value results 
either from a larger than one local p- value in the su- 
percritical regime n > 1 or from a small (smaller than 
1) local p- value before the exponential growth regime. 

Such a small p-value precursor was first proposed 
empirically by Liu [1984], who studied several after- 
shocks sequences of moderate earthquakes that have 
been followed by a large earthquake. He proposed 
that a p- value smaller than 1 is a signature of a fore- 
shock sequence, whereas p > 1 is associated with nor- 
mal aftershock sequences with a single mainshock in 
the past. He suggested that p- values close to one char- 
acterize double-mainshock sequences. These empiri- 
cal rules are part of the earthquake prediction method 
used in China [Liu, 1984; Zhang et al, 1999]. The 
small precursory p-value has been used with other 
precursors to predict the occurrence of a M — 6.4 
earthquake in China following an other M = 6.4 
earthquake three months later [Zhang et al., 1999]. 
A precursor associated with a small global p-value 
has also been observed by Bowman [1997] for a se- 
quence in Australia. In 1987, several M = 4 — 5 
earthquakes occurred in a region that was not seis- 
mically active before, and triggered a large number 
of aftershocks characterized by an abnormally low p- 
value of 0.3. A sequence of three M > 6 occurred one 
year later, followed by an aftershock sequence with a 
more standard p- value of 1.1. Simeonova and Solakov 
[1999] have also reported a very low p- value of 0.5, for 
one sequence of aftershocks in Bulgaria, that was fol- 
lowed one year latter by a larger earthquake. The first 
part of the aftershock sequence was well fitted by a 
modified Omori law, and then a significant deviation 
occurred with an abnormally high aftershock rate by 
comparison with the prior trend. This departure from 
an Omori law is expected from our results for an af- 
tershock sequence in the super-critical regime and the 
very low value of the exponent p can be interpreted as 
the apparent exponent within the cross-over from the 
l/t 1 ^ 6 decay (14) at early times to the exponential 
explosion (17) at times t > t* (see Figure 3). 

In addition to the small precursory p-value pre- 
dicted in the regime n > 1, we have shown that this 
regime is also characterized by a decrease of the ap- 
parent global p-value with time. Such a decrease of 
p-value has also been identified as a precursor by Liu 
[1984]. 

Other patterns may be a signature of the super- 
critical regime. The relative precursory quiescence 
suggested by Drakatos [2000] may also be explained 
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by our results. In contrast to the "absolute" quies- 
cence which detects changes in the background seis- 
micity after removing the aftershocks from the catalog 
[e.g. Wyss and Habermann, 1988], the "relative" qui- 
escence [Matsu'ura, 1986; Drakatos, 2000] takes into 
account the aftershocks and detects changes in seismic 
activity after a large mainshock by comparison with 
the usual Omori law decay of aftershocks. Drakatos 
[2000] studied several aftershock sequences in Greece 
which contains large aftershocks, i.e. aftershock with 
magnitude no smaller than M — 1.2, where M is the 
mainshock magnitude. For each sequence, he fitted 
the aftershock sequence by a modified Omori law up 
to the time of the large aftershock using a maximum 
likelihood method. He found that large aftershocks 
were often preceded by a relative quiescence by com- 
parison with an Omori law, with an increase of the 
seismicity rate just before the large mainshock oc- 
currence. Such a departure from an Omori law is 
predicted by our results in the super-critical regime. 
Indeed, in the super-critical regime, large aftershocks 
are likely to occur when the earthquake rate N(t) 
changes from an Omori law to the exponential explo- 
sion for times close to t* . 

To illustrate this concept, we have performed a 
simulation of the ETAS model in the super-critical 
regime and have applied the same procedure as used 
by Drakatos [2000] to fit the synthetic aftershock se- 
quence by an Omori law up to the time of the first 
large aftershock. The parameters of the synthetic cat- 
alog are K = 0.024, m = 0, c = 0.001, a = 0.5, 
b = 0.8 and 9 = 0.2, yielding n = 1.27 and t* = 4.6. 
Figure 6 represents the cumulative aftershock number 
as a function of time for the synthetic catalog and the 
fit with a modified Omori law. From this figure, we 
see a clear relative seismic quiescence, as defined by a 
cumulative aftershock number smaller than that pre- 
dicted by the fit. The aftershock activity recovers the 
level predicted by the fit at the time of the large after- 
shock. All theses results are similar to those obtained 
by Drakatos [2000]. 

In the case n > 1, our results predict an expo- 
nential increase of the seismicity rate at large times. 
Because we assume that the magnitude distribution 
is independent of time, the same exponential accel- 
eration is expected for both the cumulative energy 
release and the cumulative number of earthquakes. 
Sykes and Jaume [1990] found that several large 
earthquakes in the San- Francisco Bay area where pre- 
ceded by an acceleration of the cumulative energy re- 
lease that can be fitted by an exponential function, as 



predicted by our results. In laboratory experiments 
of rupture, several studies have also observed an ex- 
ponential acceleration of the seismic energy release 
before the macroscopic rupture [Scholz, 1968; Mered- 
ith et al, 1990; Main et at, 1992]. 

More recently, many studies have reported an ac- 
celeration of seismicity prior to great events (see 
[Sammis and Sornette, 2001; Vere-Jones et al., 2001] 
for reviews) but they used a power-law instead of an 
exponential law to fit the acceleration of seismicity. 
A power-law increase of the seismicity before rupture 
is predicted by several statistical models of rupture in 
heterogeneous media, which consider the global rup- 
ture or the great earthquake as a critical point (see 
[Sornette, 2000a] for a review). Note that it is of- 
ten difficult to distinguish in real data an exponential 
increase from a power-law increase, especially with a 
small number of points and for times far from the rup- 
ture time. No systematic study has been undertaken 
that compares these two laws to test if the accelera- 
tion of the seismicity is better fitted by a power-law 
rather than by an exponential law (see however [Jo- 
hansen et al, 1996]). 

We have stressed that the ETAS model is funda- 
mentally a mean field approximation (branching pro- 
cess) which neglects "loops", i.e., multiple interac- 
tions (see Figure 1). An important consequence of 
this approximation is that the super-critical regime 
cannot lead to a growth rate faster than exponential. 
Indeed, recall that an exponential growth is charac- 
terized by a time derivative of the number of events 
proportional to the number of events dN / dt — N/t* , 
i.e., is fundamental a linear process. In a sequel to 
the present work [Sornette and Helmstetter, 2001], 
we show however that for b < a, the impact of the 
largest earthquake induces an effective nonlincarity 
which leads to a faster-than-exponential growth rate, 
possibly leading to a finite-time singularity [Sam- 
mis and Sornette, 2002]. A faster-than-exponential 
growth rate may also be obtained by introducing mul- 
tiple interactions between earthquakes and positive 
feedback: rather than the linear law dN/dt — N/t* 
expressing the condition that each "daughter" has 
only one "mother" , we may expect an effective law 
dN/dt ~ N s , with 5 > 1 providing a measure of the 
effective number of ancestors impacting directly on 
the birth of a daughter. We may thus expect that an 
improvement of the ETAS model beyond the "mean- 
field" approximation would lead to power law accel- 
eration of seismicity. 

Other precursory patterns may also be related to 
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the super-critical regime: they comprise the precur- 
sory earthquake swarm or burst of aftershocks [Evi- 
son, 1977; Keilis-Borok et al, 1980a, 1980b; Molchan 
et al., 1990; Evison and Rhoades, 1999]. Swarms are 
earthquake sequences characterized by high cluster- 
ing in space and time and the occurrence of several 
large events with magnitude larger than M — 1, where 
M is the magnitude of the largest event. A burst of 
aftershocks is a sequence of one or more mainshocks 
with abnormally large number of aftershocks at the 
beginning of their aftershock sequences [Keilis-Borok 
et at, 1980a]. From our results, an abnormally high 
aftershock rate or a sequence with several large events 
are expected in the super-critical regime. 

Temporal change of n-value and transition 
from one regime to the other one 

It is often reported that the b and p values vary 
in space and time [e.g., Smith, 1981; Guo and Ogata, 
1995, 1997; Wiemer and Katsumata, 1999]. We have 
documented that a part of the observed variation of 
the exponent p may not be genuine but result from an 
inadequate parameterization of a more complex real- 
ity. Because n and t* are function of b, p and the other 
ETAS parameters, we expect the fundamental param- 
eters of the ETAS model, namely n and t* , to vary 
significantly in space and time. The branching ratio 
n plays the role of a "control" parameter quantifying 
the distance from the critical point n = 1 between the 
sub-critical and the super-critical regime; t* is a cross- 
over time and is sensitive to details of the systems. As 
a consequence, it is very reasonable to expect that the 
Earth's crust will change from the sub-critical to the 
super-critical regime and vice-versa, as a function of 
time and location. 

Equation (3) shows that the branching ratio n is 
a decreasing function of b. Accordingly, this may ra- 
tionalize the observation that large earthquakes are 
sometimes preceded by a decrease of the 6-value [e.g. 
Smith, 1981]. A decrease of the b- value leads to an 
increase of the n-value, that can move the seismic- 
ity from the sub-critical to the super-critical regime, 
and thus increase the probability to observe a large 
earthquake. Other ETAS parameters (a, K, p and 
c) may also change in time and move the seismic- 
ity from one regime to the other one. Ogata [1989] 
measured the ETAS model parameters before and af- 
ter the 1984 Western Nagano Prefecture earthquake 
(M = 6.8). He found that the seismic activity preced- 
ing the mainshock was characterized by a lower b, c, 
K parameters and local p values than the seismicity 



following the mainshock. He also obtained a larger 
a-value for the seismicity preceding the mainshock. 
All these changes of parameters, except the change in 
K, lead to a larger n-value before the mainshock than 
after. Before the mainshock, n is in principle infinite 
because the local p-value is smaller than one. As we 
already discussed, this corresponds to an explosive 
super-critical regime of growing seismicity. After the 
mainshock, we find n = 0.92 and t* = 10 6 days, using 
the determination of the ETAS parameters. The seis- 
micity has thus changed from a super-critical regime 
before the mainshock to a sub-critical regime after the 
mainshock. 

Conclusion 

We have provided analytical solutions of the ETAS 
model, which describes foreshocks, aftershocks and 
mainshocks on the same footing. Each event trig- 
gers an aftershock sequence with a rate that decays 
according to the local Omori law with an exponent 
p = 1 + 6. The number of aftershocks per event 
increases with its magnitude. We suggest that the 
Earth's crust at a given time and location may be 
characterized by its branching ratio n, quantifying its 
regime. We propose that n is a fundamental parame- 
ter for understanding and characterizing the organiza- 
tion of the seismicity within the Earth's crust. In the 
sub-critical regime (n < 1), the global rate of after- 
shocks (including secondary aftershocks) decays with 
the time from the mainshock with a decay law dif- 
ferent from the local Omori law. We find a crossover 
from an Omori exponent 1 — 9 for t < t* to 1 + 9 
for t > t* . The modified Omori law is thus only an 
approximation of the decay law of the global after- 
shock sequence. In the super-critical regime (n > 1 
and > 0), we find a novel transition from an Omori 
decay law with an exponent 1 — 9 at early times to an 
explosive exponential increase of the seismicity rate 
at large times. The case 9 < leads to an infinite 
n-value, due to the slow decay with time of the local 
Omori law. In this case, we find a transition from an 
Omori law with exponent 1 — \9\ similar to the local 
law, to an exponential increase at large times, with a 
crossover time r different from the characteristic time 
t* found in the case 9 > 0. These results can rational- 
ize many of the stylized facts reported for foreshock 
and aftershock sequences, such as the suggestion that 
a small p- value may be a precursor of a large earth- 
quake, the relative seismic quiescence preceding large 
aftershocks, the positive correlation between b and 
p- values, the observation that great earthquakes are 
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sometimes preceded by a decrease of 6-value and the 
acceleration of the seismicity preceding great earth- 
quakes. 

Finally, we would like to mention that our analysis 
can be generalized to various other choices of the local 
Omori law and of the magnitude distribution. The 
ETAS model can also be extended to describe the 
spatial distribution of the seismicity [Helmstetter and 
Sornette, 2002]. 
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Appendix: technical derivation of the 
analytical solution 

In this appendix, we provide the technical derivation of 
the results used in the main text for the sub-critical and 
super- critical regimes. We start from equation (9). 

General derivation for 9 > 

The integral over r is the convolution of X m i with 4> m ' ■ 
Since there is an origin of time and we have a convolution 
operator, the natural tool is the Laplace transform f(f3) = 
J +o ° /(i)e _/3t dt. Applying the Laplace transform to (9) 
yields 

X m W) = S(f3, m) +P(m) / dm' 4> m , {(3) X m , (/?) . ( Al) 

J mo 

where the r.h.s. has used the convolution theorem that 
the Laplace transform of a convolution of two functions is 
the product of the Laplace transform of the two functions. 
Let us now apply the integral operator dm 4> m (P) on 
both sides of (Al) and define 



AOS) 



f 

•J mi 



dm 4> m {/3) A m (/3) 



p oc 

Q(fi) = dm 4> m (f3) P(m) , 

and 

f oo 

S{(3) = dm 4> m (J3) S((3,m) 

Then, expression (Al) yields 

A(/3)=«S(/3)+Q(/3)A(/3) , 
whose solution is 

m 



A(/3) 



1 - Q(/3) 



(A2) 
(A3) 

(A4) 
(A5) 
(A6) 



This expression gives X m (t) after inversion of the integral 
operator J m dm (j> m (f3) and of the Laplace transform. 

The key quantity controlling the dependence of X m (t) 

is 



dm I0"("o) P(r 



dt- 



-0Ct 



(A7) 

obtained by replacing the expression of 4>m(t) defined in 

(1) and normalizing t/c -> t. Using P(m) = ln(10) b iQ-K^-m a ) 

we obtain 

Q{f3) = n R{/3c), (A8) 
where we have used the expression (3) of n and defined 



where 
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1) 



l+fl 



T(a,x) = 



= 9 e p 13" r(-6», 13) = l-e' 3 13" r(i-6 
(A9) 



(A10) 



is the (complementary) incomplete Gamma function [Abramowitz 
and Stegun, 1964] and we have used r(l+a, x) = aF(a, x)+ 
x a e~ x obtained by integration by part. Using the expan- 
sion of the incomplete Gamma function [Olver, 1974] 



r(o,x) = r(o)-^ 



k a+/e 



(-i) fc * 

k\{a + k) ' 



for a > , (All) 



we obtain 



R(f3) = l-T(l-0) /3 6 +Y~Z~§ P+O{I3 1+0 ,I3 2 , f3 2+6 ,(3 3 , ...) . 

(A12) 

It is possible, using the full expansion of the incomplete 
Gamma function, to estimate the value of A(/3) when 
the second term — ^ (3 of the expansion cannot be ne- 
glected anymore compared with the term proportional 
to f3 e . Thus, the expansion (A12) using the first two 
terms only R{/3) = 1 — T(l — 6) (3 6 becomes invalid for 
(3 > [r(l - 0)(1 - 6)] 1/{1 - e \ i.e., for times smaller than 
[r(2-6»)]- 1/(1 - e) . For all practical purpose, this is a small 
value and we can use safely the expansion (A12) in the fol- 
lowing calculations. 

Let us now make explicit A(/3): 



k r 

J m, 



dm 10' 



a(m — m a ) 
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Jo 



dt Xm(t) e 



-pt 



(A13) 

Using the definition of X(t) given by (11) and the fac- 
torization of the times and magnitudes in (A13), we obtain 



where 



X((3) = nR{(3c)X{f3) , 

f-OC 

X(f3) = dt X(t) e-^ 
Jo 



(A14) 
(A15) 
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Replacing (A14) in (A6) gives 



A(/?) = 



508) 



nR(/3c) (1 - nR(/3c)) 



(A16) 



When a great earthquake occurs at the origin of time 
t = with magnitude M, S(t,rn) = S(t) 8{m - M), ex- 
pression (A4) gives 

S (f ] ) = Je 10 Q(M " mo) i?(/3c) . (A17) 
Thus, expression (A16) becomes 

b- a io«( M - m o) 



A(/3) 



& {l-nR(/3c)) 



(A18) 



The dependence of A(/3) on /3 is uniquely controlled by 
the denominator 1 — nR(fJc). 

The sub-critical regime n < 1 

The analysis proceeds exactly as in [Sornette and Sor- 
nette, 1999]. For < 9 < 1, and for small /3 (large times), 
X(J3) given by (A18) is 



A(/J) = 



So 



So 



1 - n[l - d(/3c)«] (1 - n) U + (/%*)<> 



(A19) 

where i* is defined by (12) and the external source term So 
is defined by (13). We retrieve equation (13) of [Sornette 
and Sornette, 1999] with the correspondence to — * c. 
Two cases must be distinguished. 

• (3t* < 1 corresponds to t > t* by identifying as usual 

the dual variable f) to t in the Laplace transform with 1/t. 

i _ 



In this case, we can expand , , , which leads to 



So 
1-ra 



[1 - 09t*) fl ]- 



(A20) 



We recognize the Laplace transform of a power law of 
exponent 9, i.e. 



"(*)' 



So 



r(fl)(i- 



for t > t* 



(A21) 



• For t < t* , /3t* > 1 and (A19) can be written with a 
good approximation as 



At<*.(/?) 



So 



(1 - n)(pt*Y 



0- 



(A22) 



In other words, t* is indeed the transition time at which 
the "short-time" regime \t<t» (t) crosses over to the "long- 
time" regime A t>t *(i). 

We now calculate the full expression of A(t) valid at all 
times. We expand 

1 11 1 °° 

*)« + i = T&Ff (8t*)-o + i = TfWy ) 



-ke 



Thus, by taking the inverse Laplace transform 
S„ 1 



(A24) 



A(t) = 



1 — n 2m 



e& t j2(-i) k (pt*y( k+i) 
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(A25) 



The inverse Laplace transform of p-^ 1 ^ is t (fc+1) ^ 1 /r((fc+ 
1)9). This allows us to write 



k (t/n k 



r((k + i)e) 



(A26) 



Expression (A26) provides the solution that describes the 
cross-over from the Omori's law (A23) at early 

times to the l/t 1+e Omori's law (A21) at large times. The 



r((fc+i)fl) 



is a series representation of the 



series £)1 (-1 
special Fox function [Glockle and Nonnenmacher, 1993] 
and it is also related to the generalized Mittag-Lefner func- 
tion. 

For large times t >> t* , a direct numerical evalua- 
tion of A(i) from equation (A26) is impossible due to the 
very slow convergence of the series. The pade summation 
method [Bender and Orzag, 1978] can be used to improve 
the convergence of this series and to evaluate numerically 
(A26) for all times. 

The super-critical regime n > 1 

We can analyze this regime by putting n > 1 in (A18) 
which can be written under a form similar to (A19): 



A(/?) = 



So 



So 



Denoting F(z) = / () + °° dt e _< we see that J+ 00 dt e -73 ' t^ 1 
T{z)f3~ z . Comparing with (A22), we thus get 



A t <t.(t) 



So 



t*- 



r(0)(l-n) t 1 - 6 



for t < t* 



(A23) 



We verify the self-consistency of the two solutions 

At>t* (t) and \t<t» (t) by checking that A t >t* (t*) = \t<t» (t*). 



(l-nR((3c)) dn((3c) e - (n - 1) (n - 1) \(/3t*) e 

(A27) 

In the second and third equalities of (A27), we have used 
the small /3-expansion (A12) of R({3c) valid for < 6 < 1. 

At early times c < t < t* , i.e., fit* > 1, X(/3) w 
(„-i)(/3t*) e wm °h is the Laplace transform of (A23): thus, 
the early decay rate of aftershocks is the same ~ 1/t 1-8 
as for the sub-critical regime (A23). However, as time in- 
creases, the dual /3 of t decreases and A(/3) grows faster 
than ~ (/3c) _ due to the presence of the negative term 
— (n — 1). This can be seen as an apparent exponent 
#a PP > increasing progressively such that dn(f3t*) — 1 w 
C(f3t*) 6ap v, where C is a constant. Note that 6> app > 9 for 
the pure power law C(/3c) 9app to mimic the acceleration 
induced by the negative correction — (n— 1). The seismic 
rate will thus decay approximately as ~ i/t 1 - e a PP (*)_ 
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The large time behavior is controlled by the pole at 
(3 = \/t* of \(/3). Close to l/t*, 



AOS): 



So 



(n - 1)6» fit* - 1 
The inverse Laplace transform is thus 



(A28) 



A(i) = (2m) 



c- 



d/3 e' 3 * X(/3) 
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(n- l)t*0 



(A29) 

exhibiting the exponential growth at large times. Expres- 
sion (12) shows that ~ |1 — n|*. Thus, as expected, 
the exponential growth disappears as n — > 1 + . 

We now calculate the full expression of A(i) valid at all 
times. We expand 
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Thus 



A(t) 



(n -1) 2m . ' ^ yt ' 
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The inverse Laplace transform of l//3 (fc+1)e is t (fc+1)£ '- 1 / r ((fe+ 
1)6>). This allows us to write 



A(t) = 



So t* 



(i/t*) fce 



r((fc + i)e) 



(A32) 



Expression (A32) provides the solution that describes the 
cross-over from the Omori's law at early times to 

the exponential growth at large times. Note that the solu- 
tion (A32) can be obtained directly from (A26) by remov- 
ing the alternating sign (— l) k in the sum. The solution 
(A32) retrieves the two regimes discussed before. 

1. For t < t*, the sum in (A32) is close to 1/T(9), 
which leads to 



X(t) 



So t* 



r(0)(n-l) t 1 - 



(A33) 



2. For t > t* , the sum dominates. The sum is very 
similar to the series expansion of e t ^ t and is actu- 
ally proportional to e*/* for large t. This result is 
obvious for 9 = 1 since the series expansion becomes 
identical to that of e (,/t . This can be justifies for 
other values of 9 as follows. For 9 — > 0, the dis- 
crete sum transforms into a continuous integral of 
the type 



f 

Jo 



dx t x /T(x) 



(A34) 



A saddle-node approximation, performed using the 
Stirling approximation (which already gives a very 



good precision for small z) T(z) £s \J2-k e z z z ? , 
shows that the saddle-node of the integrant occurs 



t/t* 



For 



for x ~ t/t*, which then gives A(t) ~ e 
arbitrary 9, we can use the Poisson's summation 
rule 

+°° r + oo +°° r+oo 

f( r ) = / du f(u) + / du f{u) cos[27rqit] 

J — CO i J — OO 

r— — oo q— 1 

(A35) 

on the function defined by 



V(r9 + ( 



for r > 



(A36) 



and f(r) = for r < 0. The left-hand-side of (A35) 
is nothing but the semi-infinite sum in (A32). The 
first term in the right-hand-side retrieves the inte- 
gral (A34) encountered for the case 9 — > 0. This 
term thus contributes a term proportional to e'/* . 
All the other terms contribute negative powers of t 
and are thus negligible compared to the exponen- 
tial for t > t*. This can be seen from the fact that 
each term with q > 1 is similar to the sum in (A26) 
for the subcritical case with alternating signs. The 
larger q is, the faster is the frequency of The leading 
dependence X(t) ~ e J ^* valid for any < 9 < 1 re- 
trieves the limiting behavior already given in (A29) 
from a different approach for large times t » t*. 
It has also been proved rigorously in [Ramselaar, 
1990]. 

Case 9 < corresponding to a local Omori's 
law exponent p < 1 

The general equation (9) still holds in this case and 
the general derivation starting with (Al) up to (A6) still 
applies. The key quantity controlling the dependence of 
\ m (t) is still Q(/3) defined by (A7). Writing 9 = -|0|, we 
have 

Q{f3) = n R'iPc), (A37) 
where no is defined by (4) and 



R'((3) 



dt 



(t + l) 1 -!*! 



,/3) (A38) 



where F(a, x) is the (complementary) incomplete Gamma 
function defined by (A10). Using the exact expansion 
(All), we obtain 



Q(P)=noe f3c (Pc)-^ r(|0|)-^ 



\e\+k 



(~l) fc (fie) 
k\ (\9\+k) 



\ k = 

(A39) 

For small /3's (i.e., large times), expression (A39) has the 
following leading behavior 

Q(P) = no Y(\9\) (/3 C )-I fl ' - j^J + n Q T(\9\) (/3c) 1 -"" +h.o.t. 

(A40) 
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where h.o.t. stands for higher-order terms in the expansion 
in increasing powers of /3c. 

The source term S(f3) in the denominator of A(/3) given 
by (A6) is now given by 

S{0) = K c w i Q ( M - m °> R'(f3c) . (A41) 

Expression (A6) for A(/3) then yields 

So 



A(/3) 



(A42) 



1 - Q(f3c) ' 

where R'(/3c) is given by (A38), no is defined by (4) and 
So is defined by (13). The dependence of A(/3) on /3 is 
uniquely controlled by the denominator 1 — Q{/3c) = 1 — 
n () R'(f3c). 

Using (A40), we get the leading behavior for small /3c 

So 



AGS) = 



So 1 



(A43) 



(l + ^)(l-(/3r)-|e|) 

where the characteristic time r is given by (19). 
At early times c<t < r, (/3r) _|e| < 1 so that 

m^TT^r, {l + (rpr W ) ■ (A44) 
+ \e\ > 

By applying the inverse Laplace transform, the constant 
term contributes a Dirac function S(t) which is irrelevant 
as the calculation is valid only for t > c. The other term 
(jf3)- w gives 



A(t) = 



So 



-\»\ 



(A45) 



(l + ffi)r(|0|) 

The early time behavior of \(t) is thus similar to the local 
Omori law l/i 1 " 16 * 1 . 

Similarly to the super-critical case n > 1 of the regime 
9 > 0, the long time dependence of the regime 6 < is 
controlled by a simple pole ft* — i . 

Thus, the long-time seismicity is given by 

So 



A(t) 



(A46) 



We can also calculate the full expression of \(i) valid 
at all times t > c. We expand 



1 



1-(/3t)-I»I 



X> t >" 



fc|9| 



(A47) 



Removing the constant term, which by the inverse Laplace 
transform contributes a Dirac function S(t) which is irrel- 
evant as the calculation is valid only for t c, we get 



A(t) 



So 



d/3 e' 3 ' ^ (pry 

k=l 



k\e\ 
(A48) 



The inverse Laplace transform of l/(3 kW is t^- 1 /T(k\0\). 
This allows us to write 



S 1 00 

m = tttW) i E 



r(fc|fl|) 



(A49) 



Expression (A49) provides the solution that describes the 
cross-over from the local Omori law l/i 1- ' 9 ' at early times 
to the exponential growth at large times. 
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Figure 1. Schematic representation of the branching process associated with the ETAS model defined by (1) and 
(9). In this example, the thickest dashed line is the time arrow associated with the main shock indicated as T'. 
This main shock triggered five aftershocks denoted '11', '12', '13', '14' and '15' whose magnitudes are proportional 
to the length of their vertical lines (their position above or below the thickest dashed line is arbitrary and chosen 
to ensure a better visibility of the diagram). The aftershock '11' triggered three aftershocks denoted '111', '112' 
and '113'. The aftershock '12' triggered four aftershocks denoted '121', '122', '123' and '124'. The aftershock '13' 
triggered a single aftershock denoted '131'. The aftershock '14' also triggered a single aftershock denoted '141'. 
The aftershock '15' did not trigger any aftershock. The observable catalog is the superposition of all these events 
which are projected on the thick dashed line at the bottom of the figure, keeping the thickness as a code for the 
generation number of each event. 
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Figure 2. Seismicity rate N(t) in the sub-critical regime with n = 0.95. The noisy black line represents the 
seismicity rate obtained for a synthetic catalog generated using K = 0.024, M = 6.8, mo = 0, c = 0.001 day, 
a = 0.5, b = 1.0 and 9 = 0.2, giving the characteristic time is t* — 4500 days. The local Omori law with exponent 
p = 1 + 9 = 1.2 is shown for reference (dotted line). The analytical solution (16) is shown as the thick line. The 
two dashed lines represents the asymptotic solutions (14) for t < t* and (15) for t > t*. 
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Figure 3. Seismicity rate N(t) in the super-critical regime. Same legend as in Figure 2. The synthetic catalog 
was generated using the same parameters as for Figure 2, except for a lowest b- value of b = 0.75 and a smallest 
mainshock magnitude M = 6, leading to a branching number n = 1.43 and a characteristic time t* = 0.85 day. 
The analytical solution (thick line) is calculated from equation (18). The two dashed lines correspond to the 
approximative analytical solutions (14) and (17). 
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Figure 4. Seismicity rate N(t) in the case < corresponding to a local Omori's law exponent p < 1. Same legend 
as in Figure 2. The synthetic catalog was generated using K = 0.02, M = 7, too = 0, c = 0.01 day a — 0.5, b = 1.0 
and 9 = —0.1, giving the characteristic time is r = 10 5 days. The analytical solution (thick line) is calculated from 
equation (22). The two dashed lines correspond to the approximative analytical solutions (20) and (21). 
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Figure 5. Comparison between the three decay laws of aftershock sequences: Omori law with p = 0.7 (dashed 
line), stretched exponential with q = 1.3 and to = 10 days (thin black line) and our analytical solution in the 
sub-critical regime (16) for 9 = q = 1 — p = 0.3 and t* = to = 10 days (solid gray line). At early times t << t*, the 
three functions are similar and decay as t~° . At large times, the stretched exponential function and the analytical 
solution of the ETAS model decay more rapidly that the Omori law. For times up to t = 10 t* , the stretched 
exponential function is a good approximation of the ETAS model solution, and describes the transition from a 
power law decay at early times to a faster decay law. 
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Figure 6. Cumulative aftershock number in the super-critical regime from a synthetic catalog generated using 
branching ratio n = 1.27, 9 = 0.2 and t* = 4.6 days. The mainshock magnitude is M = 7.0. The thin line is 
a fit by an Omori law evaluated for time before the occurrence of the first M > 6.0 aftershock. This fit gives 
an apparent global p-value of 0.58. Relative seismic quiescence (by comparison with an Omori law) is observed 
before the occurrence of the M — 6.0 aftershock, due to the transition from an Omori law decay with exponent 
p = 1 — 9 = 0.8 for time t « t* to an exponential increase of the seismicity rate for time t >> t* . 
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Table 1. ETAS parameters, branching ratio n and characteristic time t* for the sequences studied by Ogata [1989, 
1992]. We have computed n and t* using equations (3) and (12) from the ETAS parameters K, a, c, p = 1 + and /i 
calculated by Ogata [1989, 1992] using a maximum likelihood method. For most sequences, we have assumed b = 1 to 
evaluate n and t* because 6-value is not given in [Ogata, 1989, 1992]. Thus, there is a large uncertainty in the n and t* 
values in the case where a is close to 1. 



Reference seismicity data M b [i K c p a n t* 

day -1 day day 



Ogata 


[1989] 


Japan, 1895-1980 


6.0 


1 


Ogata 


[1989] 


Rat-Island 1963-1982 


4.7 


1 


Ogata 


[1989] 


Nagano, 1978-1986 


2.5 





Ogata 


[1989] 


Nagano aftershocks, 1986 


2.9 


1 


Ogata 


[1992] 


worldwide shallow earthquakes 


7.0 


1 


Ogata 


[1992] 


Central Aleutian, 10 years 


4.7 


1 


Ogata 


[1992] 


Tohoku, 95 years 


6.0 


1 


Ogata 


[1992] 


Tokachi-Oki aftershocks, 1 year 


4.8 


1 


Ogata 


[1992] 


Niigata aftershocks, 150 days 


4.0 


1 


Ogata 


[1992] 


Niigata aftershocks, 150 days 


2.5 


1 


Ogata 


[1992] 


Izu Islands, 55 years 


4.0 


1 


Ogata 


[1992] 


Izu Peninsula, 7 years 


2.5 


1 


Ogata 


[1992] 


Off cast cost of Izu, 33 days 


2.9 


1 


Ogata 


[1992] 


Matsushiro swarm, 20 years 


3.9 


1 


Ogata 


[1992] 


Kanto, 1904-1916 


5.4 


1 


Ogata 


[1992] 


Kanto, 1916-1923 


5.4 


1 


Ogata 


[1992] 


Hachijo, 1938-1969 


5.4 


1 


Ogata 


[1992] 


Hachijo, 1969-1973 


5.4 


1 


Ogata 


[1992] 


Tonankai, 1933-1939 


5.2 


1 


Ogata 


[1992] 


Tonankai, 1939-1944 


5.2 


1 


Ogata 


[1992] 


Tokachi, 1926-1945 


5.0 


1 


Ogata 


[1992] 


Tokachi, 1945-1952 


5.0 


1 


Ogata 


[1992] 


Tokachi, 1952-1961 


5.0 


1 


Ogata 


[1992] 


Tokachi, 1961-1968 


5.0 


1 






0.005 


0.087 


0.02 


1.0 


0.7 


Inf 


a 





0.0 


0.072 


0.167 


1.35 


0.63 


1.04 


4600 


9 


0.021 


0.008 


0.017 


0.85 


0.94 


Inf 


b 


2 


0.0 


0.032 


0.038 


1.14 


0.73 


0.92 


4.10 6 





0.019 


0.018 


0.21 


1.03 


0.53 


1.49 


10 17 





0.008 


0.042 


0.03 


1.13 


0.62 


1.34 


2200 





0.0054 


0.98 


0.02 


1.0 


0.70 


Inf 


a 





0.14 


0.015 


0.23 


1.28 


0.98 


4.03 


1.5 





0.075 


0.0005 


0.15 


1.37 


1.26 


Inf 


b 





0.47 


0.0002 


1.10 


1.72 


1.34 


Inf 


b 





0.0038 


0.062 


0.012 


1.143 


0.155 


0.96 


10 8 





0.022 


0.035 


0.003 


1.35 


0.17 


0.91 


7.3 





0.59 


0.016 


0.009 


1.73 


0.31 


1.00 


346. 





0.0006 


0.092 


0.13 


1.14 


0.27 


1.21 


2200 





0.028 


0.010 


0.010 


1.00 


0.62 


Inf 


a 





0.025 


0.001 


0.010 


1.02 


1.31 


Inf 


b 





0.013 


0.008 


0.004 


1.02 


0.85 


3.0 


5.10 6 





0.016 


0.001 


0.013 


1.00 


1.11 


Inf 


a 





0.050 


0.010 


0.065 


1.02 


0.90 


5.28 


4.10 3 





0.031 


0.009 


0.011 


1.01 


0.83 


5.54 


10 7 





0.047 


0.013 


0.065 


1.32 


0.83 


0.57 


0.40 





0.041 


5.20 


11.6 


3.50 


1.37 


Inf 


b 





0.032 


0.021 


0.059 


1.10 


0.72 


0.99 


10 22 





0.014 


0.014 


0.005 


0.86 


0.43 


Inf 


7.10 5 



a t* cannot be evaluated because p = 1 
h t* cannot be evaluated because a < b 
c t is given instead of t* because 9 < 



